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This  paper  deals  with  the  Large-Eddy  Simulation  of  a  cryogenic  flame  issued  from  a 
LOX-CH4  shear  coaxial  injector.  The  operating  pressure  is  above  the  critical  pressure  for 
both  propellants  but  oxygen  is  injected  below  its  critical  temperature.  Such  a  state  is 
referred  to  as  transcritical  and  is  representative  of  the  extreme  conditions  which  prevail 
in  liquid  rocket  engines.  Transcritical  flows  also  exhibits  large  variations  of  density  and 
strong  departure  from  the  perfect  gas  behaviour.  To  handle  this  problem,  the  solver  uses 
a  hybrid  upwind-central  scheme  able  to  capture  the  sharp  density  gradients  and  the  real 
gas  thermodynamics  is  modeled  using  a  cubic  equation  of  state.  Finite-rate  chemistry  is 
modeled  using  a  modified  one  step  CH4/air  mechanism.  Then,  the  filtered  reaction  rate 
for  each  species  is  closed  by  correcting  the  resolved  reaction  rate  with  a  subgrid  turbulent 
contribution.  Two  simple  forms  of  this  closure  have  been  tested  and  compared.  Qualitative 
comparisons  with  experimental  data  show  that  the  LES  is  able  to  capture  the  turbulent 
structure  of  the  flame.  A  dominant  diffusion  mode  of  combustion  is  reported  and  the  flame 
is  observed  to  be  attached  to  the  LOX  post. 


I.  Introduction 

Large  Eddy  Simulation  (LES)  is  becoming  a  standard  tool  to  study  turbulent  combustion.  The  application 
of  this  method  to  model  flames  in  liquid  rocket  engines  remains  however  a  challenging  task.  The  difficulty  of 
such  a  simulations  is  mainly  due  the  transcritical  conditions  which  usually  prevail  in  the  system:  propellants 
are  injected  at  cryogenic  temperature  below  the  critical  temperature  (T  <  Tc)  while  pressure  exceeds  the 
critical  value  (P  >  Pc).  Under  this  regime,  surface  tension  and  latent  heat  become  null.* 1 * *  Atomization  and 
vaporization  are  replaced  by  continuous  mixing  processes  dominated  by  diffusion  and  turbulent  convection.2-4 
Moreover  gas  and  liquid  can  not  be  distinguished:  transcritical  fluids  may  be  very  dense  while  keeping  gas 
like-properties.5 

During  the  last  two  decades,  several  experimental  studies  have  been  carried  out  on  cryogenic  flames6-9 
and  they  have  shown  important  differences  between  subcritical  and  supercritical  combustion.  At  subcritical 
pressure,  the  competition  between  dynamic  forces  and  surface  tension  entails  the  liquid  core  break-up  and 
generates  a  spray.  Then  atomization  is  followed  by  vaporization  and  combustion.  For  this  regime,  the  flame 
stabilization  mechanism  depends  of  the  fuel  (CH4  or  H2)  and  of  the  injection  conditions.  The  flame  may 
be  lifted  in  the  Oxidizer /Fuel  shear  layer  or  it  may  be  anchored  at  the  LOX  post.9  Above  the  supercritical 
pressure,  droplets  are  replaced  by  small  finger-like  structures  which  dissolve  in  the  light  medium  surrounding 
the  transcritical  jet.7  Under  these  conditions,  the  flame  attachment  on  the  injector  tip  is  systematically 
observed. 

Substantial  progress  have  also  been  obtained  in  the  modeling  and  the  simulation  of  transcritical  com¬ 
bustion.  Harstad  et  al.10  proposed  a  computationally  efficient  method  based  on  the  departure  function 
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and  the  Peng- Robinson  Equation  of  state  to  describe  thermodynamic  properties  of  fuel/oxidant  mixtures. 
This  method  has  been  widely  applied  in  CFD  and  has  shown  good  accuracy.  Large  Eddy  Simulation  and 
Direct  numerical  simulations,  mainly  carried  out  on  mixing  layers  and  coaxial  configurations  for  LOX/H2 
flames,4,11,12  have  extended  the  understanding  of  the  processes  which  control  transcritical  combustion.  By 
computing  a  sector  of  the  single  element  experiment  of  Oschwald  et  al,7  Oefelein13  was  able  to  provide  a 
detailed  characterization  of  the  LOX/H2  flame  near  field.  He  also  identified  the  primary  holding  mechanism 
leading  to  the  flame  attachment  on  the  injector  tip. 

This  paper  proposes  a  LES  strategy  to  solve  transcritical  oxy- combust  ion  of  methane.  The  choice  of 
methane  is  driven  by  the  recent  interest  for  this  fuel  as  an  alternative  to  hydrogen.  Even  if  it  has  not  been 
widely  utilized  for  spacecraft  applications,  many  arguments  suggest  that  the  next  generations  of  rockets  will 
include  engines  operating  with  cryogenic  methane.  Its  diffusivity  is  lower  and  its  liquefaction  temperature 
is  higher  than  hydrogen.  This  would  allow  to  design  lighter  tanks  and  by  extension  might  increase  the  pay 
load  of  the  future  space  vehicles. 

The  combustion  of  methane  at  high  pressure  also  raises  new  questions  and  new  challenges.  First,  the 
density  and  the  critical  temperature  of  this  gas  are  high  compared  to  hydrogen.  This  allows  a  large  range 
of  injection  condition  and  by  extension  a  large  range  of  combustion  regimes  that  CFD  tools  must  be  able 
to  handle.  The  experiment  carried  out  on  a  shear  coaxial  injector  by  Singla  et  al.14  illustrates  this  point. 
They  have  shown  that  when  methane  and  oxygen  are  both  injected  below  the  critical  temperature,  the  flame 
exhibits  a  large  expansion  angle.  Moreover,  combustion  takes  place  in  two  distinct  regions:  in  the  inner 
and  the  outer  shear  layers  of  the  coaxial  jet.  The  same  experience  repeated  with  supercritical  methane 
provides  one  zone  of  combustion  in  the  shear  layer  between  the  oxygen  and  the  fuel  streams  and  a  smaller 
angle  of  expansion  is  observed.  Finally,  if  detailed  mechanisms  can  be  now  applied  to  model  finite  rate 
chemistry  in  LES  of  hydrogen- oxygen  flames,15  their  equivalents  for  methane  remain  too  expensive  to  be 
utilized  with  high  fidelity  simulation  technics.  For  comparison,  the  detailed  high  pressure  hydrogen-oxygen 
mechanism  proposed  by  Conaire  et  al16  consists  of  21  elementary  steps  among  8  species  while  the  GRI-Mech 
3.017  mechanism  contains  325  steps  and  53  species.  This  issue  adds  complexity  to  the  open  question  of 
the  the  reaction  rate  closure.  The  flamelet  approach18  which  assumes  an  infinitely  fast  chemistry  is  a  quite 
inexpensive  method  which  can  circumvents  this  problem.  But  for  many  rocket  applications  where  the  flame 
may  experience  high  strain  rate  an  local  extinction,  this  method  becomes  inappropriate  since  the  chemical 
time  scales  and  the  flow  time  scales  may  be  comparable.  Finite  rate  chemistry  closures  such  as  the  Linear 
Eddy  Model  (LEM)19  would  be  more  suitable  to  simulate  this  kind  of  configuration. 

In  the  present  study,  a  Large  Eddy  Simulation  of  the  LOX/CH4  Mascotte  test  rig  (Version  4) 14  has  been 
carried  out  for  the  operating  point  G2.  Recent  numerical  studies  of  this  case  using  RANS20  and  LES21  with 
infinitely  fast  combustion  models  have  provided  a  good  representation  of  the  experimental  flame  and  have 
shown  strong  real  gas  effects  on  its  structure.  This  paper  first  presents  the  formulation  applied  to  model 
the  LES  equations,  the  thermodynamics  and  the  transport  properties.  A  simple  LES  finite  rate  chemistry 
approach  is  also  proposed  to  close  the  combustion  reaction  rate.  Then,  the  experimental  configuration 
and  the  numerical  setup  are  detailed.  Finally,  the  calculation  of  the  transcritical  flame  is  discussed  and 
qualitatively  compared  with  the  experiment. 

II.  Formulation 


II. A.  LES  closure 


The  current  solver  deals  with  a  Favre-filtered  version  of  the  compressible,  multi-species  and  unsteady  Navier- 
Stokes  equations  in  conservative  form.  Following  Erlebacher  et  al.,22  the  flow  variables  are  decomposed  by  a 
spatial  filtering  operation  (denoted  by  an  overbar  /)  and  two  separate  fields  are  obtained:  the  unresolved,  or 
sub-grid,  scale  and  the  resolved,  or  super-grid,  scale  represented  by  a  tilde  (/  =  ^=~).  The  filtered  equations, 
respectively  for  mass  conservation  (1),  momentum  conservation  (2),  energy  conservation  (3)  and  species 
conservation  (4),  read: 


dp  dpUi 
dt  dxi 


(1) 


dpUi 

dt 


d 

dxj 


(pUiUj  +  pSij  -  Tij  +  r-fs) 


=  0 


(2) 
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(3) 


~^t~  +  -&x-  (peTUi  +  pUi  +  *Aik  “  uiTii  +  HiSS  +  <T*SSS)  =  0 
^jr  +  (pYkUi  +  Ji,k  +  Y-,t  +  0^)  =  dik  for  k  =  1  •  •  •  Ns  (4) 


The  closure  of  the  system  above  is  an  issue.  Ns  +  5  equations  have  been  written  to  obtain  the  Ns  +  5 
conservative  variables, but  several  other  quantities  have  been  introduced:  the  filtered  pressure  p,  the  heat 
flux  Qi  IK  (in  its  “Irving-Kirkwood”  form  including  the  enthalpy  flux  by  mass  diffusion),  the  mass  flux 
and  several  sub-grid  terms  (r2sgs ,  i72sgs,  cqgs,  Y^sgs  and  6^ ) .  The  subgrid  heat  flux  g\9S  and  the  subgrid 
species  diffusive  flux  are  usually  neglected  because  they  have  been  shown  to  have  a  small  magnitude  at 
high  Reynolds  number.23,24 

The  mass  and  heat  diffusion  fluxes  are  solely  based  on  Fick’s  and  Fourier’s  types  of  diffusion: 


- dT 


Ns 


Ns 


Qi, ik  —  ~ ^  ^  hkYkVi,k  +  Qll 


sgs 

k 


k=l 


k= 1 


Ji,k  —  pYkVi,k  —  P'^k,m 


(5) 

(6) 


with  the  approximation  of  evaluating  the  diffusion  coefficient  of  each  species  into  the  mixture  instead  of 
considering  each  individual  binary  diffusions.  Thus  the  cross-diffusion  Dufour  and  Soret  terms  are  neglected 
in  this  study.  Their  significance  under  super-critical  conditions  is  not  always  demonstrated25  and  their 
inclusion  in  a  LES  formulation  would  introduce  additional  unclosed  terms  that  have  not  been  studied  yet. 

Like  cqgs  and  #®gs,  the  sub-grid  enthalpy  flux  due  to  mass  diffusion  Q ^  is  neglected  as  a  first  approxi¬ 
mation.  The  sub-grid  enthalpy  flux  iLzsgs,  the  sub-grid  stress  tensor  r2sgs  and  the  sub-grid  species  flux  Y^sgs 
are  closed  using  a  gradient  diffusion  approach  and  a  sub-grid  eddy  viscosity  that  uses  the  turbulent  kinetic 
energy. 

The  momentum  closure  will  define  the  eddy  viscosity.  An  additional  transport  equation  is  needed  for  this 
turbulent  kinetic  energy  ksgs  =  5  (tijtii  —  UiUi).  Neglecting  compressibility  effects  and  assuming  a  simple 
gradient  diffusion  model  for  the  sub-grid  transport,  this  equation  reads: 


dpk3gs 

dt 


+  —  (pkagsUi)  =  Psgs  -  esgs 

OXi 


_  vt  dksgs 
PPrt  dxi 


(7) 


cscra  £) —  \/(kS&S}  3 

The  expressions  for  the  production  Psgs  =  —  and  dissipation  esgs  =  Cep v  ^ —  terms  leave  only 

two  model  coefficients,  Cu  and  Ce  to  complete  the  sub-grid  closure.  In  this  study,  Cv  and  Ce  are  constant, 
respectively  equal  to  0.067  and  0.916.  These  values  were  evaluated  in  earlier  studies27  without  any  a  priori 
assumptions  of  ideal  gas  behavior  but  assuming  incompressible  flows.  The  turbulent  Prandtl  and  Schmidt 
numbers  are  also  assumed  constant  with,  respectively,  values  of  0.9  and  0.7. 


II. B.  Thermodynamics  and  transport  properties 

Among  the  cubic  equation  of  state  listed  by  Poling  et  al.,1  the  Peng- Robinson  equation  of  state  is  the  one 
chosen  for  this  study: 

=  RuT _ A 

P  Vm-B  V*+  2 VmB  -  P2  1  ’ 

where  Vm  is  the  molar  volume,  Ru  is  the  universal  gas  constant.  The  first  term  models  the  repulsive 

force  that  molecule  exert  on  each  other  at  short  distance.  The  parameter  B  is  proportional  to  the  actual 
volume  of  the  molecule.  The  second  term  y2  +2^  b-b2  m°dels  the  long  range  attractive  forces  between  the 
molecules  such  as  electrostatic  forces,  polarization  or  London  dispersion  forces. 

The  transport  properties  are  also  affected  by  the  transcritical  conditions.  It  is  required  to  apply  robust 
models  which  are  able  to  handle  large  variation  of  the  fluid  properties  from  a  dense  cryogenic  fluid  in  the 
LOX  flow  to  a  light  gas  in  the  flame.  For  this  study,  a  correction  of  the  Chung  model  for  high  pressure  has 
been  applied:  the  viscosity  fi  and  the  thermal  conductivity  k  are  computed  using  the  Corresponding  States 
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Method.  The  basic  viscosity  equation  for  a  rigid  non-attracting  sphere28  are  solved,  including  the  correction 
by  Chung  et  al.  for  the  shape  and  polarity.29  The  diffusion  coefficients  are  modeled  with  a  correction  on 
the  basic  equation  proposed  by  Fuller  et  al.30 

II. C.  Reduced  Mechanism 

If  detailed  mechanisms  for  CH4  combustion  such  as  GRI-Mech  3.0  are  now  available,  they  involve  too  many 
species  and  reaction  steps  to  be  applicable  in  LES.  Only  a  reduced  mechanism  with  a  few  steps  and  species 
can  be  reasonably  used  for  this  type  of  simulation.  However  they  are  usually  built  for  Air/CH4  flames  and 
they  match  a  small  range  of  thermodynamic  conditions  (pressure,  initial  temperature...).  For  this  study, 
CH4  has  to  burn  with  pure  oxygen  under  high  pressure.  Therefore,  correction  of  the  one  step  mechanism 
of  Westbrook  and  Dryer31  is  introduced  here  to  predict  adequate  adiabatic  flame  temperature  and  flame 
speed  under  these  conditions.  The  original  mechanism  is  composed  of  Nr  =  1  irreversible  reaction  involving 
Ns  =  4  species  CH4,  02,  C02,  H20: 


CH4  +  202  ->  C02  +  2 H20 


(9) 


It  will  be  called  WD1  in  the  rest  of  the  paper.  Preliminary  simulations  of  a  CH4/02  laminar  premixed 
flame  at  ip  =  1  with  an  initial  temperature  Tq  =  300 K  were  carried  out  using  Cantera.  The  pressure  was 
Po  =  5 .6MPa  which  corresponds  to  the  ambient  pressure  of  the  case  G2  (See  Table  2).  The  real  gas  equation 
of  state  was  replaced  by  the  perfect  gas  equation.  The  final  composition  in  the  burnt  gas  is  easy  to  compute 
for  WD1: 

Xcha  =  Xo2  =  0,  XCo2  =  ^  0  =  1/3,  Xh2o  =  — =  2/3  (10) 

+  2  (p  +  2 

Moreover,  the  final  adiabatic  temperature  Tad  is  given  by  applying  the  first  thermodynamic  law  to  the 
system: 

VpQ  -  X°i)Ah°f?  +  Xi  [Tad  cmP)idT  -  X°  [T°  cmP)idT  =  0  (11) 

2=1  JTa  JTref 

The  resolution  of  Eq.  (11)  provides  Tad  =  5072AT.  Table  1  compares  the  flame  speed  and  the  adiabatic  flame 
obtained  with  WD1  and  GRI-Mech  3.0.  WD1  presents  big  discrepancies  with  the  detailed  kinetic  chemistry 
which  shows  this  mechanism  is  not  adapted  to  oxycombustion. 

A  simple  method  to  improve  WD1  is  presented  here.  The  new  reduced  mechanism  will  be  called  WDlox. 
First,  The  adiabatic  temperature  is  overpredicted  because  the  one  step  mechanism  does  not  take  account 
of  dissociations.  These  reactions  play  an  important  role  in  oxycombustion.  This  problem  can  be  solved  by 
adding  a  term  in  the  enthalpy  equation  (Eq.  (11))  to  model  the  chemical  enthalpy  of  products  which  should 
occur  during  the  reaction  (CO,  H  ,H2...). 


A  h0,rn 

f, dissociation 


-  X°i)Ah°f  f  +  X,  [Tad  cmP!idT  -  Xf  [T°  ■ 

i— i  JTq  JTref 


jdT  =  0 


(12) 


Ahj  dissociation  can  a^so  expressed  in  function  of  the  chemical  enthalpy  variation  using  a  coefficient  o:,/ISS : 

4  fTad  rT0 

(l-adiss)J2(Xi-X°i)Ah°f7  +  Xi  cmPiidT-X°  cmP)idT  =  0  (13) 

2—i  JT0  JTref 

Eq.  (13)  can  be  solved  using  the  adiabatic  temperature  predicted  by  GRI-Mech  3.0.  It  gives  OLdiss  =  0.31.  In 
practice,  the  effect  of  dissociation  can  be  added  by  multiplying  the  enthalpy  of  formation  for  each  species  by 
1  —  &diss-  In  Cantera  and  LESLIE,  thermochemistry  is  computed  using  the  NASA  polynomial  interpolation: 


Hi/RT  =  oi  +  a2T/  2  +  a3T2/ 3  +  a4T3/4  +  a5T4/ 5  +  a6/T  (14) 

According  to  Eq.  14,  the  enthalpy  of  formation  for  each  species  i  is  equal  to  clq.  In  order  to  obtain  the  correct 
adiabatic  temperature,  this  coefficient  is  multiplied  by  1  —  QLdiss-  The  other  quantities,  crnP:i  and  Si  are  not 
affected  by  this  modification  since  they  don’t  depend  of  clq. 
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Then  the  flame  speed  sl  must  also  be  changed.  The  asymptotic  analysis  of  Zeldovich,  Frank-Kamenetski 
and  Von  Karman  (ZFK)32  demonstrates  that  sl  varies  like  the  square  root  of  the  reaction  rate  coefficient 
Af.  Accordingly,  a  way  to  match  the  correct  flame  speed  consist  in  modifying  Af. 

Af correction  =  Af  f - ^  (15) 

\SLWD1oxJ 

where  slgri  is  the  flame  speed  obtained  with  GRI-Mech  3.0  and  slwdiox  is  the  flame  speed  given  by 
WDlox  without  any  correction  of  Af.  The  third  line  of  Table  1  presents  the  adiabatic  temperatures  and 
the  flame  speed  obtained  with  WDlox.  The  corrected  mechanism  predicts  these  quantities  with  less  than 
5%  error.  However,  WDlox  presents  some  limitations.  Figure  1  shows  that  the  agreement  between  WDlox 
and  GRI-Mech  3.0  is  not  perfect.  The  heat  release  maximal  value  is  too  high  while  it  is  underestimated  for 
temperatures  lower  than  2000K.  The  ignition  temperature  is  around  1500K  whereas  GRI-Mech  3.0  predicts 
combustion  up  to  750K.  This  issue  can  create  extinctions  when  hot  reactants  are  mixed  with  some  colder 
gas  by  turbulent  eddies. 


Tad  (K) 

sL  (m/s) 

GRI-Mech  3.0 

3584 

2.317 

WDl 

5051 

11.11 

WDlox 

3755 

2.283 

Table  1.  Adiabatic  temperature  and  flame  speed  obtained  with  Cantera 


Figure  1.  Premixed  flame  heat  release  versus  temperature:  comparison  of  WDlox  with  GRI-Mech  3.0 

II. D.  LES  Combustion  Model 

The  last  term  to  close  in  Eq.  4  is  the  LES  filtered  reaction  rate  4//..  The  strategy  chosen  here  consists  in 
splitting  ujk  between  a  large  scale  and  a  subgrid  contribution: 

Zjk  =  Fuk  {p,  T,  Y?)  +  Gusk9S  (16) 

The  large  scale  reaction  rate  (i>k  (^p,  T,  Y))  is  directly  expressed  from  the  Arrhenius  laws  of  the  chemical 

mechanism  applied  to  the  resolved  field.  The  subgrid  reaction  rate  ujsk9s  depends  of  the  subgrid  scale  pa¬ 
rameters,  such  as  the  subgrid  kinetic  energy,  and  must  be  modeled.  The  terms  F  and  G  are  the  blending 
functions  which  depend  of  the  local  flow  conditions.  Two  simple  forms  of  these  functions  have  been  tested: 
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LAM  F  =  1  and  G  =  0:  the  subgrid  term  is  neglected.  Mixing  is  assumed  perfect  at  the  subgrid  level. 


MAX  F  =  H(\uk\  >  \(^l9S\)  and  G  =  1  —  F:  where  H  is  the  Heaviside  function.  This  model  is  equivalent 
to  pick  the  maximum  value  between  the  laminar  reaction  rate  and  the  subgrid  reaction  rate.  The 
objective  of  this  model  is  to  replace  the  large  scale  reaction  rate  by  the  usk9S  only  in  the  zones  where 
the  turbulent  effects  are  high  and  where  the  flame  can  be  diluted  by  cold  gas  which  would  entail  the 
underestimation  of  the  reaction  rate  by  WDlox  and  in  some  situations  would  generate  extinctions. 
Therefore  the  model  MAX  forces  combustion  in  locations  where  oxygen  and  methane  are  in  contact 
but  where  the  reaction  rate  predicted  by  the  Arrhenius  law  is  too  low. 


The  procedure  to  compute  the  subgrid  reaction  rate  is  based  on  the  Eddy  Dissipation  Combustion  model 
(EDC).  This  approach  is  an  evolution  of  the  Eddy  Break-Up  model  (EBU)33,34  and  was  originally  developed 
by  Magnussen  and  Hjertager 35-37  for  RANS.  For  each  reaction,  a  subgrid  mixing  reaction  rate  is  computed: 


•  If  the  reaction  kinetics  goes  forward  (formation  of  products)  then  the  mixing  is  limited  by  the  deficient 
reactant: 


-mm 


k 


1  sgs 


wk 


(17) 


If  the  reaction  kinetics  goes  backward  (formation  of  reactants)  then  the  mixing  is  limited  by  the 
deficient  product: 

P  •  (  Yk  \ 


-mm 


1 sgs 


.MT. 


(18) 


where  ukr.od  and  z/£|ac  are  the  stoichiometric  coefficients  and  rsgs  is  the  subgrid  turbulent  time  scale  which 
is  estimated  from  the  subgrid  kinetic  energy  Ksgs  and  the  cell  volume  Vceii- 


Tsgs 


1/3 


cell 


VK, 


sgs 


(19) 


Finally  the  subgrid  reaction  rate  is  expressed  for  each  species  k: 


UJ 


sgs 


nreac 

=  wk  ]T  (i 


fprod 

'kj 


reac\  • 

ukj  ) 


J=1 


(20) 


This  formulation  is  very  similar  to  the  LES  subgrid  model  proposed  by  Fureby  and  Moller.23  However  the 
interpretation  of  Cosks  and  Cok  T,  Y^j  is  different.  Fureby  an  Moller  picked  the  minimum  value  of  each 

contribution  considering  that  the  slowest  phenomenon  (mixing  or  chemistry)  controls  the  reaction  rate.  This 
approach  is  consistent  for  RANS  or  coarse  LES  where  a  non  negligible  part  of  the  turbulent  mixing  occurs 
at  the  subgrid  scale.  However  for  fine  grids,  it  would  asymptotically  yields  to  a  reaction  rate  equal  to  zero 
(' Tsgs  =  0).  When  most  of  the  turbulent  mixing  scales  are  resolved  by  LES,  the  reaction  rate  is  actually 

controlled  by  the  chemistry:  | uk  \  »  \Cj^s\.  Therefore  the  decomposition  of  Eq.  16  is  more 

accurate  and  converges  to  the  correct  DNS  flame. 


III.  Experimental  and  numerical  setup 
III. A.  The  MASCOTTE  V04  test  case 

This  paper  focuses  on  the  simulation  of  the  LOX/CH4  single-element  shear  coaxial  injector  investigated 
by  Singla  et  al.14  on  the  Mascotte  test  rig.  This  small  scale  thrust  chamber  is  operated  by  ONERA.  A 
detailed  description  of  the  experimental  setup  and  of  the  instrumentation  is  given  in  previous  papers.14,38 
The  combustor  is  composed  of  a  square  chamber  equipped  with  visualization  windows  on  its  four  sides.  The 
mean  flame  structure  was  obtained  from  Abel  transform  of  OH  and  CH  LIF  emission  images.  Therefore  only 
qualitative  comparison  between  experiment  and  LES  will  be  possible.  Four  operating  points  were  studied, 
but  only  the  case  G2  at  pressure  above  the  critical  point  can  be  considered  in  this  study.  The  flow  conditions 
for  this  configuration  are  listed  in  Table  2.  Methane  is  injected  as  a  supercritical  gas  and  the  oxygen  remains 
transcritical. 
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TO  (g/s) 

T(  K) 

P  (MPa) 

7 rr 

02 

44.4 

85 

5.61 

1.11 

CH4 

143.1 

288 

5.61 

1.24 

Table  2.  Experimental  operating  conditions  for  the  case  G2.  The  reduce  pressure  7rr  designates  the  ratio  between  the 
operating  pressure  P  and  the  critical  pressure  Pc 


III.B.  Numerical  implementation 

Even  if  surface  tension  and  latent  heat  of  vaporization  vanish,  transcritical  flows  still  present  very  large 
density  gradient  due  to  mass  and  heat  transfer  from  dense  to  light  fluids. A  previous  study  on  L0X/H2 
flames39,40  have  pointed  the  limitations  of  central  scheme  to  predict  such  large  density  gradients  and  have 
shown  the  necessity  of  running  a  hybrid  central-upwind  scheme  instead  of  a  pure  central  scheme. 

Consequently,  the  current  fully  compressible  flow  solver  uses  a  finite-volume  scheme  with  a  second- 
order  time-accurate  predictor-corrector  integration  and  a  second-order  accurate  hybrid  solver  for  spatial 
integration.  This  finite  volume  hybrid  solver  alternates  between  a  second-order  central  scheme  and  a  third- 
order  accurate  MUSCL,  an  upwind-biased  scheme.  A  dynamic  and  local  switch,41  based  on  pressure  and 
density  gradients,  determines  (at  each  time  step  and  computational  face)  which  scheme  to  use.  The  MUSCL 
reconstruction  technique  is  used  alongside  an  approximate  Riemann  solver  (specifically,  Harten-Lax-Van 
Leer  (HLL)  from  Genin  and  Menon,42  with  HLL  contact /solver  modifications  by  Toro  et  ah, 43  as  well  as 
the  monotonized  central  limiter,  to  enforce  the  total  variation  diminishing  condition.  This  hybrid  scheme 
allows  the  capture  of  the  large  density  gradients  typically  found  near  the  injection  plane  while  keeping  the 
required  grid  resolution  reasonable  and  keeping  the  less  dissipative  central  scheme  in  the  far  field  in  order 
to  accurately  model  the  turbulence.  Characteristic  non-reflecting  boundary  conditions  are  employed  at  fuel 
and  oxidizer  inlets  and  at  the  computational  domain  outlet.  No-slip  adiabatic  conditions  are  imposed  at 
walls. 

III. C.  Computational  grids 

The  computational  domain  is  shown  on  Fig.  2a.  It  corresponds  to  the  first  150  mm  of  the  full  experimental 
combustion  chamber  and  includes  the  final  part  of  the  single  shear  coaxial  injector.  The  reduction  of 
the  combustion  is  motivated  by  the  very  small  length  of  the  flame  which  is  lower  than  50  mm  for  all 
operating  points.  The  length  of  the  computational  domain  is  sufficient  to  avoid  any  interaction  between  the 
reacting  zone  and  the  outlet.  Finally,  this  operation  enables  to  improve  the  resolution  without  increasing  the 
computation  cost.  The  simulations  have  been  performed  on  a  multiblock  grid.  The  outer  cylindrical  block 
contains  765  x  233  x  129  grid  points,  while  the  inner  butterfly  Cartesian  block  contains  765  x  33  x  33  grid 
points  (Fig.  2c)  The  finest  resolution  is  located  in  the  injector  posttip  region  ( Ax  =  30  p m).  The  resolved 
kinetic  energy  spectrum  in  this  region,  shown  in  Figure  3,  demonstrates  the  recovery  of  the  Kolmogorov 
—5/3  spectrum,  indicating  sufficient  resolution  for  a  proper  LES  model. 

IV.  Results 

This  section  discusses  the  results  obtained  for  the  two  LES  models  described  in  section  II. D.  At  this  stage, 
only  the  LES  of  the  model  MAX  has  been  completed.  The  laminar  chemistry  model  LAM  experienced  strong 
extinction  issues  and  could  not  be  run  enough  to  provide  converged  statistics.  By  default,  most  of  the  results 
presented  in  this  papers  will  therefore  refer  to  the  subgrid  closure  MAX.  Explicit  references  to  LAM  will  be 
carried  out  when  the  two  models  are  compared. 

IV.  A.  Flame  and  flow  Structure 

To  illustrate  the  turbulent  flow  structure,  Figure  4a  shows  iso-contours  of  Q  criterion  colored  with  the 
velocity  magnitude.  The  density  isosurface  p  =  150  kg/s  is  also  presented  in  purple  to  visualize  the  dense 
core.  Annular  coherent  eddies  occur  in  the  shear  layer  between  the  methane  stream  and  the  quiescent 
environment.  They  entail  gas  from  the  low  velocity  zone  to  the  high  speed  jet  .  These  structures  wrinkle 
the  transcritical  oxygen  core  and  break  down  to  generate  a  highly  turbulent  flow  .  The  LES  temperature 
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(a)  Overall  view 


(b)  Injector  near-field  (c)  Butterfly  grid 


Figure  2.  Views  of  the  grid  for  the  Mascotte  configuration 
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Figure  3.  Resolved  velocity  and  kinetic  energy  spectrum  in  the  injector  near  field. 


isosurface  T  =  1700LC  is  compared  with  a  visualization  of  the  experimental  flame  on  Figure  4b.  This 
qualitative  comparison  suggests  a  good  agreement  between  LES  and  experiment.  The  flame  is  anchored  on 
the  tip  between  the  CH4  and  the  LOX  injection.  The  reaction  zone  is  a  short  cone  (around  6  cm)  starting 
with  an  expansion  angle  lower  than  10  °  and  abruptly  terminating  with  an  angle  of  20°. 

Figure  5  presents  the  average  (top)  and  the  RMS  (bottom)  axial  velocity  field.  The  white  line  designates 
the  temperature  isocontour  T  =  1300  K  and  helps  to  visualize  the  combustion  zone.  The  flame  strongly 
delays  the  transition  from  a  coaxial  jet  to  a  single  jet.  Instead  of  converging  towards  the  center  line  of  the 
chamber,  the  methane  stream  is  deflected  towards  the  combustor  walls  owing  to  the  gas  expansion  generated 
by  the  flame.  A  small  zone  of  negative  axial  velocity  is  also  visible  in  the  vicinity  of  the  LOX  post.  This 
wake  plays  an  important  role  in  the  holding  mechanism  of  the  flame.13  This  creates  recirculation  of  the  hot 
products  and  stabilizes  the  flame  by  promoting  the  vaporization  and  the  mixing  of  the  dense  LOX  core. 

This  phenomenon  is  illustrated  on  Figures  6  and  7,  which  show  instantaneous  snapshots  of  the  flow  in 
the  vicinity  of  the  injector.  For  each  subgrid  model,  a  strong  backflow  of  hot  products  anchors  the  flame  on 
the  tip  wall.  But  some  differences  are  noticeable.  The  laminar  model  presents  a  preferred  flame  stabilization 
on  the  oxidizer  side  of  the  injector  while  it  is  on  the  fuel  side  for  the  model  MAX.  But  the  direct  anchoring 
of  the  flame  is  actually  unphysical.  OH-PLIF  measurements  carried  out  by  Singla  et  al38  have  shown  that 
the  flame  is  lifted  at  a  very  short  distance  of  the  wall.  Catching  this  phenomenon  with  LES  is  beyond  the 
scope  of  our  study.  The  adiabatic  wall  should  be  replaced  by  an  isothermal  boundary  condition.  This  would 
also  require  a  much  finer  resolution  in  the  vicinity  of  the  injector  and  a  detailed  mechanism  able  to  predict 
low  temperature  kinetics. 

IV.B.  Combustion  regime 

Figure  8a  presents  scatter  plots  of  species  mass  fraction  (CH4,  02  and  H20)  versus  the  mixture  fraction  z: 

sYCH4  ~  Yq2  +  Yg 2 

z  sy0  i  yo 

bICH4  ^  1 02 
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(a)  (b) 


Figure  4.  Flame  and  flow  structure:  (a)Isosurfaces  of  Qcriterion  (=5e8  s-2)  colored  by  the  velocity  magnitude  and 
density  isosurface  p  =  150  (purple),  (b)  top:  visualization  of  the  experimental  flame,14  bottom:  LES  temperature 
isosurface  T  =  1700 K 


Mean  axial  velocity 

10.000  40.000 


Figure  5.  Average  flow  field:  average  axial  velocity  (top),  RMS  axial  velocity  (bottom).  The  white  line  corresponds  to 
the  mean  temperature  isocontour  T  =  1300 
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(a)  XY  plane:  Temperature  field  on  the  sleeve  wall  (b)  YZ  plane:  Temperature  and  velocity  fields  in  the 

vicinity  of  the  tip  wall.  Isocontours  of  mixture  fraction 
(black  z  =  0.1,  grey  z  =  0.2,  black  z  =  0.3) 


Figure  6.  Flame  holding  mechanism  for  the  laminar  chemistry  model  LAM. 


(a)  XY  plane:  Temperature  field  on  the  sleeve  wall  (b)  YZ  plane:  Temperature  and  velocity  fields  in  the 

vicinity  of  the  tip  wall.  Isocontours  of  mixture  fraction 
(black  z  —  0.1,  grey  z  =  0.2,  black  z  =  0.3) 


Figure  7.  Flame  holding  mechanism  for  the  LES  model  MAX 
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where  s  =  4  is  the  stoichiometric  mass  ratio  and  YqH4  =  Yq2  =  1.  The  points  gather  around  straight  lines 
intersecting  at  the  stoichiometric  mixture  fraction,  z  =  0.2,  which  suggests  that  an  infinitely  fast  chemistry 
flamelet  model  is  valid  for  the  case  G2.  The  diffusion  regime  is  confirmed  by  the  flame  criterion44  on  Figure 
8b.  This  quantity  is  defined  as: 

Tf  =  YYo2.VYCh4  (22) 

Negative  values  of  T f  correspond  to  a  diffusion  flame  and  positive  values  define  zones  of  partially-premixed 
combustion.  The  criterion  shows  very  high  negative  values  along  the  stoichiometric  mixture  fraction  while 
nearly  no  positive  value  is  observed  elsewhere  in  the  flame.  This  kind  of  situation  is  specific  of  the  G2 
case  and  is  not  really  representative  of  the  usual  rocket  operating  conditions.  The  Reynolds  number  in  real 
rocket  engines  is  much  higher  and  the  flame  experiences  a  higher  amount  of  strain.  Therefore,  the  Damkohler 
number  might  take  a  small  value  and  finite  rate  chemistry  must  be  taken  into  account. 


0.1  0.2  0.3  0.4  0.5  0  6  0.7  06  0.0  1 

Mixture  Fraction  z 


(a) 


(b) 


Figure  8.  Flame  regime:  (a)  species  versus  mixture  fraction  scatter  plot,  (b)  flame  index  Ff.  1)  plane  z=0,  2)  plane 
x=0.5Dcif4,  3)  plane  x=5Dcjj4 


It  is  also  interesting  to  analyze  the  effect  of  the  models  defined  in  sections  II. D  and  II. C  on  the  LES 
temperature  field.  Figure  9  presents  two  temperature  scatter  plots  corresponding  to  the  two  subgrid  closures 
LAM  and  MAX.  The  points  are  also  colored  by  the  local  oxidizer  strain  rate  s  as  defined  by  Balakrishnan 
et  al.45  The  cloud  is  compared  with  the  simplified  Burke- Schumann  flame  structure  assuming  infinitely  fast 
chemistry  and  calorically  perfect  gas.  The  maximal  temperature  is  chosen  equal  to  the  adiabatic  temperature 
predicted  by  GRI-Mech  3.0  at  the  stoichiometry  (cf.  table  1).  For  the  two  cases,  the  evolution  of  the 
temperature  is  close  to  the  Burke- Schumann  solution  and  the  maximal  temperature  is  correctly  predicted 
which  validates  the  modified  1-step  mechanism  WDlox.  The  flame  obtained  with  the  laminar  chemistry, 
however,  exhibits  lots  of  extinctions  which  look  unphysical,  since  they  were  not  reported  by  the  experiment. 
Moreover,  by  computing  laminar  Oxygen/CH4  counterflow  flames,  Pons  et  al.46  showed  that  the  extinction 
strain  rate  at  P  =  5.6  MPa  was  sext  =  106  s-1.  Figure  10  illustrates  how  MAX  fixes  this  problem.  The 
subgrid  methane  reaction  rate  contours  are  plotted  in  color.  The  black  contours  correspond  to  values  of  the 
filtered  reaction  rate  greater  than  104  kg/m3/s.  The  subgrid  reaction  rate  is  several  order  of  magnitude  lower 
than  the  LES  reaction  rate  but  it  is  mainly  observed  on  the  periphery  of  the  flame  in  the  zones  where  the 
turbulence  start  to  strongly  wrinkle  the  flame.  The  laminar  chemistry  closure  fails  to  predict  combustion  at 
these  locations  and  yields  to  spurious  partially  premixed  zones  which  can  not  be  correctly  lighted  up  whereas 
MAX  allows  ignition  and  maintains  a  strict  diffusion  flame  regime. 
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Figure  9.  Temperature  versus  mixture  fraction  scatter  plot  colored  by  the  strain  rate  (/s):  (a)  Laminar  chemistry 

LAM,  (b)  Subgrid  closure  MAX 


Figure  10.  Methane  reaction  rate  instantaneous  field:  colored  contours  correspond  to  the  subgrid  reaction  rate 
and  black  contours  correspond  to  the  filter  reaction  rate  \ucha\  >  104  kg/m3/s. 
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IV. C.  Real  gas  effects 

The  fluid  properties  variations  across  the  flame  are  shown  in  Fig.  11.  First  Fig.  lla,b  present  the  radial 
variations  of  temperature,  density  and  species  mass  fraction  x/D=0.5.  D=0.01  m  is  the  methane  outer 
diameter.  The  density  is  divided  by  100  between  the  dense  oxygen  and  the  flame.  It  generates  very  sharp 
density  gradients  in  the  shear  layer  (0.18  <  r/D  <  0.3)  where  both  mixing  and  combustion  occur.  Then 
Fig.  lie  shows  the  compressibility  factor  Z  =  PV/RT.  This  parameter  measures  the  departure  of  the  fluid 
from  the  ideal  gas  assumption.  In  the  methane  stream  and  in  the  far  field  (r/D  >  0.6),  Z  remains  close  to 
1  while  in  the  oxygen  core  it  is  close  to  0.2.  In  this  zone,  density  predicted  by  the  Peng  Robinson  equation 
of  state  is  greater  than  1000  Kg/m3  while  the  perfect  gas  equation  of  state  would  provide  a  density  around 
230  Kg/m3.  Therefore,  even  if  the  perfect  gas  assumption  might  be  reasonable  in  most  of  the  computational 
domain,  applying  this  equation  of  state  for  this  case  would  actually  entail  important  discrepancies  on  the 
flame  expansion.  This  conclusion  is  enforced  by  Figs.  lld,e,f  which  plot  the  radial  variation  of  the  kinematic 
viscosity  and  the  Prandtl,  Lewis  and  Schmidt  numbers.  The  kinematic  viscosity  strongly  increases  in  the 
flame.  The  Prandtl  number  increases  by  a  factor  of  4  in  the  the  LOX  side  of  the  flame  to  Fuel  side  while 
the  Lewis  number  and  the  Schmidt  numbers  increases  by  a  factor  of  100.  As  observed  by  Oefelein13  for  the 
LOX/H2  flame,  the  comparison  of  these  three  numbers  shows  that  the  mass  diffusion  is  the  slowest  transport 
phenomenon  in  the  liquid  oxygen  and  is  rate  limiting  for  the  combustion  on  this  side  of  the  flame. 

V.  Conclusion 

This  paper  has  highlighted  the  results  obtained  in  the  large  eddy  simulation  of  transcritical  liquid  oxy¬ 
gen/supercritical  methane  flame.  The  solver  uses  a  hybrid  upwind  central  scheme  to  capture  the  large  density 
gradients  in  the  transcritical  flow  as  well  as  the  turbulence  created  by  the  shear  layers  and  the  break  down  of 
the  central  jet.  The  thermodynamics  of  the  liquid  oxygen  is  modeled  with  the  Peng- Robinson  cubic  equation 
of  state,  which  provides  a  good  compromise  between  cost  and  accuracy  for  Large-Eddy  Simulations.  Finite- 
rate  chemistry  is  modeled  using  a  one  step  CH4/air  which  was  modified  to  provide  realistic  temperatures 
into  the  CH4/02  flame.  Then  the  filtered  reaction  rate  for  each  species  is  closed  by  correcting  the  resolved 
reaction  rate  with  a  subgrid  turbulent  contribution.  Two  simple  forms  of  this  closure  have  been  tested.  The 
simulation  shows  a  good  qualitative  agreement  with  the  experiment.  The  two  combustion  models  achieve  to 
stabilize  the  flame  through  the  recirculation  zone  in  the  vicinity  of  the  tip  wall  but  the  model  LAM  exhibits 
local  extinctions  while  the  second  model  MAX  prevents  them.  This  simple  LES  combustion  model  is  a 
first  step  towards  more  accurate  finite  rate  approaches  able  to  handle  a  wide  range  of  combustion  regimes. 
Our  next  objective  is  to  apply  a  LEM  closure  to  this  case  with  the  same  mechanism  to  provide  a  better 
description  of  the  flame/turbulence  interaction.  A  complete  closure  using  both  LEM  and  detailed  chemistry 
remains  too  expensive  to  simulate  real  gas  combustion  in  complex  geometries.  We  are  therefore  investigating 
the  artificial  network  methodology  to  speed  up  the  chemistry  kinetics  calculation  and  circumvent  this  issue. 
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Figure  11.  Radial  variation  of  species  mass  fraction,  temperature,  density,  compressibility  factor,  kinetic  viscosity  and 
Prantdlt,  Schmidt  and  Lewis  numbers  at  an  axial  location  x/D=0.5  (D=0.01  m  is  the  methane  outer  diameter) 
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